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Abstract 

We formulate a variational fictitious-time flow which drives an initial guess torus to a torus 
invariant under given dynamics. The method is general and applies in principle to continuous 
time flows and discrete time maps in arbitrary dimension, and to both Hamiltonian and dissipative 
systems. 
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I. INTRODUCTION 



Analysis of dynamical systems in terms of invariant phase space structures provides im- 
portant insights into the behavior of physical systems. The simplest such invariants are 
equilibria, points in phase space which are stationary solutions or 0-dimensional invariants 
of the flow. They and their stable/unstable manifolds yield information about the topology 
of the flow. The role that the next class of flow invariants, periodic orbits, play in the 
topological organization of phase space and the computation of long time chaotic dynamics 
averages is well known (for an overview, see Ref. [1]). A periodic orbit is topologically a 
circle or an invariant 1-torus for a flow, and a set of discrete points or an invariant 0-torus 
for a map, embedded in a (i-dimensional phase space. Higher dimensional invariant tori also 
frequently play an important role in the dynamics; we refer the reader to Ref. [2] for some 
of the references to this literature. Invariant tori of dimension lower than the dimension of 
the dynamical flow can be normally hyperbolic [3, 1]. KAM theory implies that invariant 
tori occur in Cantor sets, and such tori play key roles in the phase-space transport [5, 6]. 
For 2-degree of freedom Hamiltonian flows (i.e., 4-dimensional phase space), 2-dimensional 
invariant tori act as barriers to diffusion through phase space, and for higher dimensional 
flows such structures act as effective barriers (Arnold diffusion). In dissipative systems such 
as Newtonian fluids, quasi-periodic motion on two or higher dimensional tori is one of the 
routes to the eventual turbulent motion [7, 8]. 

There are many methods for determining periodic orbits available in the literature [1, 9, 
10]. The lack of comparably effective methods for the determination of higher dimensional 
invariant structures has stymied the exploration of the phase spaces of high-dimensional 
flows, a focus of much recent research [2, 11]. 

Signal processing methods like frequency analysis [12, 1-3], based on the analysis of trajec- 
tories, can detect elliptic invariant tori since these tori influence signiflcantly the behavior 
of nearby trajectories. Bailout methods [14, 15] can effectively locate the elliptic regions 
in a non-integrable system, by embedding the dynamical system into a larger phase space. 
Ref. [16] describes a variational technique designed to flnd regular orbits in a phase space 
with mixed dynamics. However, these methods can only detect trajectories with non-positive 
Lyapunov exponents. They single out regular motions in a phase space but can not exactly 
determine a torus unless it is stable. Due to their relative ease of identiflcation, in partic- 
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ular cases periodic orbits are used to study invariant tori and their breakup. For example, 
in Greene's criterion approach [17, 18, 19] one studies a sequence of periodic orbits which 
converges to a given invariant torus. Such approaches have been mainly applied to the 
determination of tori of Hamiltonian systems with 2 degrees of freedom. 

Other techniques to determine invariant tori are specific to the phase space dynamics 
of the system under consideration, most often a Hamiltonian system. Early attempts like 
spectral balance method were based on the computation of quasi-periodic orbits [20, 21], 
the closure of which constitutes the invariant torus. To overcome the small divisor problems 
associated with the flow on a torus, recent research employed a geometric point of view 
and focused on the invariant torus itself. Efforts are devoted to find the solution of the so- 
called invariance condition which ensures the invariance of a parametrized object in phase 
space. Invariance conditions are functional equations for maps [22, 23, 24, 25, 26, 27] and 
first-order partial differential equations (PDEs) for flows [28, 29, 30]. These equations can 
be solved by Newton's method or Hadamard graph transform technique [1]. In view of the 
periodicity in the angle variables, Fourier transforms are widely used in the computation [31, 
32, 33, 34, 35, 36, 37]. For Hamiltonian systems, the action principle and the Hamilton- 
Jacobi equation are also frequently used in the calculation of periodic and quasi-periodic 
orbits [32, 38, 39, 40, 41, 42]. 

In this paper we introduce a method to solve the invariance condition equation and 
obtain invariant m-tori of flows and maps embedded in d-dimensional phase spaces. It 
has a variational structure which guarantees the global convergence. The method is a 
generalization of the variational "Newton descent" method originally developed to locate 
periodic orbits of flows [43, 44], which can be viewed as a variant of multi-shooting method 
in boundary value problems [45, 46, 47]. When the representative points on the guess torus 
achieve a near-continuous distribution, a PDE is derived which governs their evolution to 
a true invariant torus. In spirit, this is similar to the approach used in Ref. [32] and thus 
high accuracy is expected. However, our method is stable and thus applies to more general 
cases, including the searches for partially hyperbolic tori embedded in the chaotic regions 
of phase space. In a general dynamical system, the phase space structure can be extremely 
complex, and the global stability of our algorithm is of key importance for the the efficiency 
of our searching program. In our numerical computation, an adaptive scheme is used which 
keeps changing the step size according to the smoothness of the evolution. In addition to 
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the adaptive step size, we further speed up our searches by utihzing the continuity of the 
variational evolution equations. These points will be explained in detail in what follows. 

In Sec. II we derive the variational equation which governs the fictitious time dynamics. 
The numerical implemention of this equation is discussed in Sec. III. The method is further 
illustrated in Sec. IV through its application to the determination of 1-tori of the standard 
map, of 2-tori of a forced pendulum flow (3-dimensional phase space), of 1- and 2-tori of two 
coupled standard maps (a four dimensional symplectic map), and of 2-tori of the Kuramoto- 
Sivashinsky system (infinite dimensional phase space). In particular, we provide evidence 
that the method converges up to the threshold of existence of a given invariant torus and 
yields estimates of the critical thresholds of the breakup of invariant tori of 2-degree of 
freedom Hamiltonian systems. 

II. NEWTON DESCENT METHOD FOR INVARIANT TORI 

We start by deriving a variational fictitious time evolution equation for the determination 
of a 1-dimensional invariant torus of a d-dimensional map f : M'^ — ^> M'^. The method can be 
extended to the determination of invariant m-tori of d-dimensional maps and flows. 

A fixed point (0-dimensional invariant torus) x = f(x) is a point which is mapped into 
itself under the action of f . Likewise, a 1-dimensional invariant torus of f is a loop in M*^ 
which is mapped into itself under the action of f . If points on the invariant 1-torus are 
parametrized by a cyclic variable s G [0, 27r], with x(s) = x(s + 27r), a point x(s) is mapped 
into another point on the invariant torus 

f(x(s))=x(. + ^(.)), (1) 

where uj{s) is the local parametrization s-dependent shift. In other words, the full phase 
space dynamics f induces a 1-dimensional circle map on the invariant 1-torus 

s s + uj{s) mod27r. (2) 

We also parametrize our guess for the invariant 1-torus, the loop x(s,r), by s G [0,27r], 
with x(s, r) = x(s -|- 27r, r). Together with the "fictitious time" r, to be defined below, this 
parametrizes a continuous family of guess loops. However, for an arbitrary loop there is 
no unique definition of the shift u, as the loop is not mapped into itself under action of f. 
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Intuitively, uj should be fixed by requiring that the (i-dimensional distance vector between 
the circle map image of a point on the loop at s, and the "closest" point on the iterate of 
the loop 

F(s, r) = x(s + ^(s, r), r) - f (x(s, r)), (3) 

is minimized. For example, if the guess loop is sufficiently close to the desired invariant 
1-torus, ci;(s, r) can be fixed by intersecting the loop with a hyperplane normal to the loop 
and cutting through the image of loop f (x(s, r)). 

Compared with fixed point and periodic orbit searches for iterates of maps, the new aspect 
here is that we are searching for m-dimensional compact invariant hyper- surfaces, with 
points on such hyper- surfaces parametrized by m cyclic variables. We have encountered this 
situation already for the continuous time fiows, for which a periodic orbit p is an invariant 1- 
torus, with x(t) G p naturally parametrized by the cyclic time variable t G [0,Tp]. For other 
cyclic coordinates we are free to choose a parametrization s that best suits our purposes. 

In this exploratory foray into the world of compact higher-dimensional invariant manifolds 
we shall make the simplest choice at each turn. In particular, we are free to choose any 
parametrization s which preserves ordering of points along the invariant 1-torus, i.e. any 
circle map (2) that is strictly monotone, 1 -|- dcu/ds > 0. For an irrational rotation number 
a strictly monotone circle map can be conjugated to a constant shift, so in what follows we 
define the s parametrization dynamically, by requiring that the action of the dynamics f for 
both the guess loop and the target invariant 1-torus is rotation with a constant shift uj, 

s s + uj mod 2tt . (4) 

The invariance condition (1) with conjugate dynamics (4) has been used previously in the 
literature [33, 34]. We now design a stable scheme which yields a parametrization x(s) 
satisfying Eq. (1) together with Eq. (4). 

Following the approach of Refs. [43, 44] originally developed to locate periodic orbits of 
fiows, we now introduce the simplest cost functional that measures the average distance 
squared (3) of the guess loop from its iterate 

^^H = /^F(s,r)^ (5) 

Similar functional was used in the stochastic path extremization [48]. Here J^^[t] = J?^^[x, u] 
is a functional, as it depends on the infinity of the points x(s, r) that constitute the loop 
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for a given r. If the loop is an invariant 1-torus, = 0, otherwise JF^ > 0. At fictitious 
time r we compute cost due to the two mappings: one is the iterate f(x(s, r)) of the loop, 
and the other the circle map s i— > s + uj{t) along the loop. The fictitious time evolution 
should monotonically decrease the distance between a loop and its iterate, as measured by 
the functional JF^[t], by moving both the totality of loop points x(s,r) and modifying the 
shift uj{t). 

With constant shift circle map (4) the variation of ^^[t] under the (yet unspecified) 
fictitious time variation dr is 

^^>,^2/|(f„..).^„..)), (6) 



where 



— F(s,r) = —{s + uj{t),t) + v(s + cj(r),r) ^ ^ 



dr dr dr 

J(x(s,r))^(s,r), 



dr 



The adjustment in the loop tangent direction v is needed to redistribute points along the 
loop in order to ensure the constant shift parametrization s, and the [dxd] Jacobian matrix 
of the map J = 9f/9x moves the loop point x(s,r) in the "Newton descent" direction. 

Again we design a fictitious time flow in the space of loops by taking the simplest choice, 
in the spirit of the Newton method: 

for which jF^[x, cu] decreases exponentially with fictitious time r: 

jr2^r] = J^^[0]e~^\ (8) 

The "Newton descent" PDE (7) which evolves loop points in fictitious time r and along 
loop direction s is the main result of this paper. Written out in detail, the Newton descent 
equation for a guess loop, 

(9x , , (9x , .du,. _ , 

^{s + u,t) + —{s + u,t) — {t) (9) 

-J{x{s,t)) — {s,t) = f(x(s,r)) -x(s + u;,r). 



evolves points x(s, 0) on the r = initial guess loop to the points x(s) = x(s, oo), s t— > s + u, 
u = uj{s, oo), on the target 1-torus, provided that the r flow does not get trapped in a local 
minimum with jF^[c)o] > 0. 

The Eq. (9) can also be derived via a multi-shooting argument as has been done in 
Ref. [11]. Instead of a blind minimization of the cost functional (5), the method uses the 
vector equation (3) and its flrst derivative to flnd the zeros of the cost functional. The 
momotonicity of (8) with r ensures the global convergence. A similar argument has been 
used in the derivation of a globally convergent modifled Newton's method in Ref. [47] . 

Generalization to searches for invariant m-tori is immediate: the guess m-torus is 
parametrized by s = (si, S2, ■ ■ ■ , Sm) € [0, 27i]"^, periodic in each cyclic coordinate 



with m incommensurate shifts uj = {uji,uj2, ■ ■ ■ ,ujm) [19]- Now the fictitious time fiow (9) 
has an [d x m] invariant surface tangent tensor v. The fictitious time fiow searches (9) 
for invariant tori can also be adopted to continuous time fiows, by reducing the fiow to a 
Poincare return map on any local Poincare section which intersects the trajectories on the 
guess (m+l)-torus transversally. We will provide examples in what follows. 

In general, each tangent vector of an invariant m-torus transformation along given cyclic 
parameter Sk has a unit eigenvalue, and requires a constraint. For example, for the Jacobian 
matrix of a continuous time periodic orbit (a 1-torus) the velocity vector is an eigenvector 
with a unit eigenvalue, and Newton descent equations need to be supplemented with a 
constraint (a Poincare section) in order to determine the period of the orbit. In addition, 
if the fiow is Hamiltonian, and the invariant m-torus is located on a fixed energy surface 
H{p,(i) = E, the constraint dH/dr = is needed to ensure the conservation of the energy 
by the fictitious time dynamics. 

In case at hand, there are two alternative ways to impose the constraint: We may or may 
not fix a priori. 

(a) If we are searching for an invariant 1-torus of a fixed shift u;, the fictitious time fiow 
should not change the shift along the loop. 




for all k e 



■m 



(10) 



duj/dr = . 



(11) 
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(b) If we are searching for an invariant 1-torus of a given topology, the shift uo = uj{t) 
varies with the fictitious time r, and is to be determined simultaneously with the 1-torus 
itself. In this case we impose the phase condition [30] 

/rf.(^v(.,r)-^(s,r)^ =0, (12) 

which ensures that during the fictitious time evolution the average motion of the points along 
the loop equals zero. Empirically, for this global loop constraint the fictitious time dynamics 
is more stable than for a single-point constraint such as 5x(0, r) = 0. For m-torus, v(s, r) 
is a [d X m] tensor and Eq. (12) yields m constraints. For energy conserving Hamiltonian 
systems, one phase condition has to be replaced by the energy condition 

i- f dsWH{^{s, r)) ■ = E-^dsj H{^{s, r)) , (13) 

where a fixed E fixes the energy shell under consideration. 

The two cases are analogous to continuous time Hamiltonian flow periodic orbit con- 
straints: case (a) corresponds to fixing the period and varying the energy shell, and case (b) 
to fixing the energy and computing the period of a periodic orbit of a given topology. 

The examples of Sees. IV A, IVB and IV C illustrate the constant shift uj constraint (11); 
the examples of Fig. 4 and Sec. IV D illustrate the phase condition (12). 



III. NUMERICAL IMPLEMENTATION 

Due to the periodic boundary condition (4) it is convenient to expand the loop point x, 
the Jacobian matrix J, the map f , and the loop tangent v as a discrete Fourier series 

x(s,r) = Y^^ki^y 

k 
k 

f(x(.,r)) = Y.Mr)e''' 

k 

v(s,r) = iY.ks.k{T)e'^' (14) 

k 

(a^ = a_fc due to the reality of x(s, r), and similar relations hold for and b^), and rewrite 
the Newton descent PDE (9) as an infinite ladder of ordinary differential equations: 



Finally, the unit stability eigenvalue along the loop tangent direction v(s, r) needs to be 
eliminated by adding to (15) either the constant shift uj constraint (11), or the phase condi- 
tion (12) . in the Fourier representation the phase condition is given by Yl,k ' dsLk/dr = 

The monotone decrease with r of the functional JF^, given by (6), guarantees that the 
solution of (15) approaches a fixed point which, provided that JF^ 0, is the Fourier 
representation of the target invariant torus. 

In our numerical calculations, we represent the loop by a discrete set of points 
{x(si), ■ ■ ■ ,x(s2Ar)}. The search is initialized by a 2A^-point guess loop. The Fourier trans- 
forms of X, V and J are computed numerically, yielding M complex Fourier coefficients a^, 
bfc, and Jk, respectively. To maintain numerical accuracy, we choose M < N and set = 0, 
bfc = 0, and = for |A;| > M. We terminate the numerical integration of the fictitious 
time dynamics (18) when the distance (3) falls bellow a specified cutoff. In the Fourier 
representation, we stop when distance reaches the termination value A defined as 

max ||Ffc|| = max \bk j — ak jC^^'^l < A , (16) 

k k,j 

where akj and bkj denote the jth component of a^ and b^. 

While the algorithm is more efficient the better the initial guess, in practice it often 
works for rather inaccurate initial guesses. If the initial guess is bad, or the target invariant 
torus does not exist, the evolution diverges. Then another search is initiated, with a new 
guess. This guess torus can either be derived from the integrable limit, like the examples 
of Sees. IV A, IV B and IV C, or from numerical exploration, like the example of Sec. IV D. 
If the invariant torus is isolated or partially hyperbolic, it can be a challenging problem 
to initialize the search for an embedded invariant torus. However, once provided with a 
reasonable guess, our method is able to reliably locate the torus with relatively high accuracy. 

Another concern is related to the numerical efficiency. If we try to find a higher order 
torus (large m) in a high dimensional phase space (large d) with high accuracy, we have to 
repeatedly invert a very large [{(2M)"^d + m) x [[2M)"^d + m)] matrix when carrying out the 
integration of Eq. (15). This may constitute a major bottleneck in such calculations. In our 
numerical implementation, the matrix invertsion by the LU decomposition [47] consumes 
most of the computational time. We employ a speed-up scheme, based on the continuity 
of the evolution of Eq. (15). Once we have the LU decomposition at one step, we use it 
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to approximately invert the matrix in the next step, with accurate inversion achieved by 
iterative approximate inversion [17]. In practice, we find that one LU decomposition can 
be used for many 5t evolution steps. The more steps we go, the more iterations at each 
step are needed to get the accurate inversion. After the number of such iterations exceeds 
some fixed given maximum number, another LU decomposition is performed. The number 
of integration steps following one decomposition is an indication of the smoothness of the 
evolution, and we further accelerate our program by adjusting accordingly the step size 5t: 
the greater the number, the bigger the step size. Near the final stage of convergence, the 
evolution becomes so smooth that the step size can be brought all the way up to 5r = 1, 
recovering the full undamped Newton-Raphson step and acquiring the desired quadratic 
convergence. 

IV. EXAMPLES 

We now test the Newton descent method for determining invariant tori on a series of 
systems of increasing dimensionality: a two-dimensional area-preserving standard map, a 
Hamiltonian flow with one and half degrees of freedom (a forced pendulum), a 4-dimensional 
symplectic map (two coupled standard maps), and a dissipative PDE (the Kuramoto- 
Sivashinsky system). In the following, the representative points are uniformly distributed 
on the initial guess torus. 

A. Critical tori of the standard map 

As our first example we search for invariant 1-tori of a two-dimensional area-preserving 
map, the standard map 

Qn+i = g„ + Vn+1 mod 27r 

Pn+1 = Pn + Ksinqn , (17) 

where K is the nonlinearity parameter. For K = the map is a constant rotation in q, and 
for A" > its phase space is a mixture of KAM tori and chaotic regions. In the Fourier 
space the initial guess loop x = (g, p) and its image f (x) = [q + p + K sin q,p + K sin q) are 
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expanded as 

x(s, r) = s + ^ afc(r)e^'=' , s = (s, 0) 

k 

f(x(s,r)) = s + J2Mr)e''\ 

k 

The linear term s in Eq.(15) is needed to compensate the modulus 2tx operation on q in 
Eq.(17). Substitution into (9) yields 

/ d^k duj\ ikuj , r du ^ da.1 

= - ate"-" - Soii^ei , (18) 

where ei = (1,0). If we denote by the distance (3) on the right hand side of (18), the 
invariant torus condition for constant shift (11) is = for all k, i.e = a^e*^'^ for k 
and bo = ao + uei. 

As the first test of our variational method, we apply it to the determination of the golden- 
mean invariant torus, with shift fixed to Ug = 2it(^/5 — l)/2, and the fixed shift constraint 
(11). We use as the initial guess for the fictitious time dynamics the invariant torus of the 
linear standard map with K = and the golden-mean shift x(s, 0) = {s,ujg) , represented 
by the straight line in Fig.l. In order to test that the method works for a smooth invariant 
torus we set K = 0.5 and integrate the fictitious time dynamics (18) with 2N = 256 point 
discretization of the torus, M = 64 complex Fourier mode truncation, and A = 2 x 10~^ 
termination value (16). The resulting invariant torus is shown by the dotted line in Fig. 1. 

4.2| ^ ^ n 

4 

^3.8^ T 
3.6 

2 4 6 

q 

FIG. 1: The uj = ujg = 3.883 • • • golden mean invariant torus of the standard map (17) for K = 0.5; 
the straight line represents the initial condition. 

Next, we apply the method to a sequence of golden-mean invariant tori with increasing 
K. Numerics indicates that there exists a critical value such that when K < K^, the 
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fictitious time dynamics converges exponentially, as in (8), but for K > Kc, it diverges. The 
critical value Kc depends sensitively on the torus discretization 2N and the termination 
value A. Kc{N) computed for A = 2 x 10~^ and several values of N is 



2N 


64 


128 


256 


512 


1024 


Kc{N) 


0.34 


0.80 


0.93 


0.9656 


0.9762 



The golden-mean critical invariant torus is depicted in Fig. 2(a) for 2N = 1024 points 
discretization of the torus. Small oscillating structures in the critical torus whose resolu- 
tion would require higher frequency Fourier components are already visible. The uneven 
distribution of representative points (s parametrization's embedding into the {q,p) plane) 
along the torus indicates the drastically varying stretching rate on the invariant torus close 
to the breakup [50, 51]. Our variational method of estimating the critical Kc parameter is 
in agreement with the Greene's estimate [17] that the golden-mean invariant torus breaks 
up at the critical value Kc ~ 0.9716. Moreover, we find that for large values of 2N points 
discretization of the torus, Kc{N) approaches Kc approximately as A^~^. 



\ 

1) 246 0246 

(a) q (b) q 

FIG. 2: Invariant tori for the standard map (17) for: (a) oo = Ug at K = Kc{fy'i2) = 0.9762 close to 
the golden-mean torus critical value Kc^ termination value A = 2 x 10~^. The inset enlargement of 
the curve around q = 4.6 illustrates the fine structure of the nearly critical torus, (b) irrational shift 
uj = 2it{-k — 3) at the estimated critical value -ftrc(512) = 0.4313, termination value A = 4 x 10~^. 
2N = 1024 torus points discretization. 

As Newton descent method does not depend on the specific arithmetical properties of the 
invariant torus shift, it should work for arbitrary irrational shifts. As an example, we study 
the family of invariant tori with shift uj = 27r(7r — 3). The critical value of convergence of 
our algorithm is Kc ^ 0.4313 for 2N = 1024 and A = 4 x 10~^ The critical torus, depicted 
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on Fig. 2(b) exhibits non- uniform s-parametrization and oscillating structure, though much 
less so than the golden-mean critical torus. 

In order to assess the sensitivity of the method to the choice of the termination value 
A, we have studied its influence on the estimation of the critical Kc- For the golden-mean 
example, a decrease in the termination value to A = 10"^ for uo = uOg and 2N = 1024 points 
discretization of the torus, yields Kc = 0.6188 much smaller than the value of Kc = 0.9762 
obtained for A = 2 x 10^^. The corresponding invariant torus for A = 10~^ is depicted 
in Fig. 3(a). We notice that this torus looks much smoother than the one obtained for 
A = 2 X 10~^ (see Fig. 2(a)). Similarly, for u = 2'k{tx — 3) a decrease of the termination 
value to A = 2 X 10"^, yields a much smaller critical value K,, = 0.3004. The corresponding 
invariant torus for A = 2 x 10~^ is shown in Fig. 3(b). The points are distributed more 
evenly than in Fig. 2(b), indicating that the invariant torus obtained using this termination 
value is far from criticality. 




(a) q (b) q 

FIG. 3: The invariant tori for the standard map (17) with smaller termination values A than in 
Fig. 2, the same number of torus points 2N = 1024: (a) oj = ujg with Kc = 0.6188 and A = 10"^; 
and (b) uj = 2'k{-k - 3) with Kc = 0.3004 and A = 2 x lO"*'. 

In summary: For flxed 2N points discretization of the torus, if A is too small, then 
Kc{N) < Kc, while if A is too large, then Kc{N) > Kc- At the threshold of criticality the 
invariant torus is fractal and thus cannot be resolved by a smooth flnite Fourier truncation. 
The discrepancy between the invariant torus and its numerical discretization has a compli- 
cated influence on the flctitious time dynamics, not elucidated in this investigation. If A 
is too small, the discrepancy leads to an estimate of Kc lower than the true Kc, making 
the torus appear smoother. If A is too large, the discretization will average out the small 
features, converging to a grid beyond the critical value. With increasingly reflned 2N point- 
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discretization of the torus, the value of A needs to be chosen carefully in order to improve 
the Kc estimate. 

5| ' ' n 

4.8 '\ 

Q. / 

4.6 

4.4I ^^^r^L ^ 

2 4 6 

q 

FIG. 4: An invariant torus of the standard map (17) for K = 0.352 obtained by the fictitious time 
dynamics with the phase condition (12). The method yields shift u> 4.67857. 2N = 256 points 
discretization of the torus, termination value A = 2 x 10"^. 

So far we have determined invariant tori of the standard map by imposing a constant shift 
condition (11). An alternative is the phase condition (12) which requires that the motion of 
representative points along the torus during the fictitious time dynamics averages to zero. In 
this case the shift uj is not fixed, but is determined by the fictitious time dynamics. We test 
this condition by starting with an initial torus x(s) = (s,9ci;g/10) discretized on 2A^ = 256 
points, with termination value A = 2 x 10"^. For K = 0.352 the Newton descent method 
yields the invariant torus of the standard map shown in Fig. 4, with shift uj ^ 4.67857. 



B. A periodically forced Hamiltonian system 

As our second test case, we consider the forced pendulum 

H{p, X, t) = — e{cosx + cos(a; — t)) , (19) 

a time-dependent Hamiltonian flow with 1.5 degrees of freedom. H{p, x, t) is a periodic func- 
tion of the angle variable x and the time variable t, with dynamics on M x T^. The Poincare 
return map for the stroboscopic section t = mod 27t is a reversible area-preserving map. 
The Jacobian J required for the fictitious time dynamics (9) is evaluated by integrating 

/ 



J = AJ, A 



1 

y — e(cosx + cos(x — t)) 
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J(0) = 1 . (20) 



We apply the fixed shift condition (11) Newton descent to the determination of the invariant 
torus with the golden-mean shift uj = ujg = (y^ — l)/2. For the initial guess torus we take 
the golden-mean torus of Hamiltonian (19) with e = 0, i.e. x(s) = {s^ujg). We define ec{N) 
to be the minimum value of the parameter of the model at which the algorithm defining 
the fictitious time dynamics with 2N sampling points fails to converge at fixed A. The 
critical values edN) computed for different numbers of sampling points (termination value 
A = 2 X 10-6) are 



2N 


64 


128 


256 


512 


1024 




0.01688 


0.02312 


0.02594 


0.02750 


0.02781 



For 2N = 512 and 2N = 1024 the edN) values that we find are are close to the threshold 
£c ~ 0.02759 estimated in Ref. [52]. The invariant torus with e = 0.02781, 2N = 1024 and 
A = 2 X 10"^ shown in Fig. 5(a) exhibits non-smoothness and an uneven distribution of 
discretization points characteristic of criticality. Setting A = 10"^ leads to the invariant 
torus with the critical value estimate ic = 0.01844, displayed in Fig. 5(b). It looks smooth, 
indicating that it is far from criticality and thus that the termination value is too small. 

0.66 
0.64 

Q. 

0.62 
0.6 

2 4 6 "Id 2 4 6 

(a) X (b) X 

FIG. 5: Invariant tori of Hamiltonian (19) with to = ujg obtained by the fictitious time dynamics 
with 2N = 1024 and two different termination values: (a) A = 2 x IQ-^ yields a critical value 
Ec = 0.02781, and (b) A = lO'^ yield to an underestimate ic = 0.01844. 

C. Two coupled standard maps 

In principle, the Newton descent method is applicable to determination of invariant tori 
of arbitrary dimension for flows or maps of arbitrary dimension. In practice, one is severely 
limited by computational constraints. 
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In order to test the feasibility of the method in higher dimensions, here we consider two 
coupled standard maps [53], 

+ ei sin 9n + e-s sin(6'„, + ipn) 

On+l = On + In+l (21) 

€2 sin ipn + es sin(6'„ + ipn) 

1pn+l = + Jn+1 , 

with 4-dimensional phase space, and demonstrate that the method can determine 1- and 2- 
dimensional invariant tori. The fictitious time dynamics (15) acts on the x = {On, In, ^'n, Jn) 
phase space, with dynamics f(x) defined by (21). 

First, we apply the fixed shift (18) fictitious time dynamics to determination of the 1- 
dimensional golden mean invariant torus with shift u = Ug. For the initial guess torus we 
take the integrable case torus ei = e2 = 63 = : 

yi{s) = {s,u;g,s,ujg). (22) 

In the numerical experiment we then search for a typical 1-d invariant torus, for (arbitrarily 
chosen) small coupling values ei = 0.1, €2 = 0.15, 63 = 0.005 . 




(a) 9 (b) 9 

FIG. 6: A 1-dimensional invariant torus with shift ujg of (21) with ei = 0.1, £2 = 0.15 and 
€3 = 0.005 : {a) I — 9 projection; (b) J — 9 projection. 2N = 512 points discretization of the torus, 
termination value A = 10~^. 

The invariant torus obtained by the fictitious time dynamics in this case is shown in 
Fig. 6. Numerically 6 = ip, indicating that for this 1-dimensional torus the two phases are 
entrained. The torus appears very smooth, indicating that for the parameter values chosen 
it is far from the critical values. 
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Next, we apply the Newton descent to the determination of a 2-dimensional torus with 
non-resonant frequencies ui and uj2. In this case, we need two cychc parameters (si,S2) G 
[0, 27r]^ to locate a point on the torus. In (18) we take 





Si 




fl 


o\ 















, ei = 













1 











The initial guess is chosen as in the integrable ej = case 

S2) = (Si,UJi, S2,UJ2)- 



(23) 



In the numerical experiment we then search for (arbitrarily chosen) ei = 0.07, €2 = 0.1 and 
£3 = 0.004 2-dimensional invariant torus with (also arbitrarily chosen) frequencies Ui = ujg 
and UJ2 = vr(v^ — 1). In order to reduce the computational time, we take a rather coarse 
2N = 32 grid, with (2N)'^ = 1024 points representing the torus. 




(a) 







(b) 




FIG. 7: The 2-dimensional invariant torus of the coupled standard maps (21) with incommensurate 



frequencies wi 



and 0J2 = - 1) for ei = 0.07, 62 = 0.1 and €3 = 0.004. {2Nf = 1024 



points discretization of the torus, termination value A = 10 ^. 

Two projections of the resulting invariant torus for A = 10~^ termination value are shown 
in Fig. 7. While the 4'{si,S2) and J(si,S2) dependence on Si , S2 shown in Fig. 7 follows 
in shape the integrable case (23) dependence, the small coupling terms induce significant 
oscillations. The smoothness of the invariant torus indicates that the parameters are not 
close to the critical values. For {2N)'^ = 1024 points discretization of the torus, A can be as 
low as 5.1 X 10~^, and for {2NY = 4096, as low as 1.6 x 10~^. However, the computation 
takes at least 100 times longer, and in this exploratory study the larger (2N)'^ resolutions 
were out of reach. 
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D. Kuramoto-Sivashinsky system 



In our last example, we apply the Newton descent to determination of an invariant 2- 
torus embedded in a high- dimensional strongly contracting flow. Special tori that can be 
converted to periodic orbits in a rotating or moving frame have previously been computed 
for the complex Ginzburg-Landau equation [31], and for the 2-d PoiseuUe flow [54]. Here 
we shale determine a generic 2-torus of the Kuramoto-Sivashinsky equation [55, 56, 57] 
parametrized by the system size L, 

ut = {u^)x - u^x - u^xxx , X e[0,L]. (24) 

The Kuramoto-Sivashinsky equation describes the interfacial instabilities in a variety of 
contexts, like the flame front propagation, the two fluid model and the liquid film on an 
inclined plane. 

In the study of flame fluttering on a gas ring as the system size L increases, the "flame 
front" becomes increasingly unstable and turbulent. As shown in Refs. [8, 58], in dissipa- 
tive systems 2-dimensional tori often result from a Hopf bifurcation while 3- (or higher-) 
dimensional tori are a rare occurrence. In the following we restrict our search to the anti- 
symmetric solution space of (24) with periodic boundary conditions, i.e. u{—x,t) = —u{x,t) 
and u{x + L,t) = u{x, t), with u{x, t) Fourier-expanded as 

oo 

uix,t)= la^e'^^^, (25) 

fc=— oo 

where q = 2tx / L is the basic wavenumber and a-k = —ak G M. Accordingly, (24) becomes a 
set of ordinary differential equations : 

oo 

Ofc = {{kqf - {kqY)ak - kq Y (^mak-m ■ (26) 

m=— oo 

In the asymptotic regime of (26) for k large a^s decay faster than exponentially, so a 
finite number of a^s yields an accurate representation of the long-time dynamics. In our 
calculation, a truncation at c/ = 16 suffices for a quantitatively accurate calculation. 

In the current example, 2A^ = 128 points are used to represent the torus on the Poincare 
section ai = 0.06. Numerical experimentation indicates that for L = 40.95 trajectories 
spend significant fraction of time in a toroidal neighborhood, suggesting that a (partially 
hyperbolic?) invariant 2-torus exists at this system size: Poincare section returns of a typical 
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orbit fall close to a closed curve. The initial guess for the Newton descent is constructed by 
choosing 128 points to represent this curve and their Fourier transform is used to initialize 
the search with (18). In this case the shift uj is fixed by dynamics, and in order to compute 
it we impose the phase condition (12). 



,n -0.03 

CD -0.035 

CO 

-0.04 

' -O.O45 I 

0.08 0.1 0.12 -0.1 -0.09 -0.08 

(a) ^2 (b) S 

FIG. 8: The projections of the 2-dimensional invariant torus of (26) on the Poincare section 
ai = 0.06 with shift uj = 0.5968 for L = 40.95 : Projection on (a) (02,05) and (b) (03,00). 
The Poincare section return times are in the range T = 24.18 it 0.3. 2N = 128 torus points 
parametrization, A = 10~^ termination value. 

Fig. 8 shows two Poincare section projections, in the Fourier space, of the invariant 2- 
torus of the Kuramoto-Sivashinsky flow determined by the Newton descent method. The 
method yields the shift uj = 0.5968. Even though the invariant torus is very smooth and 
discretization points are evenly distributed, surprisingly many points are required to resolve 
the torus. For attempts with fewer discretization points, for example 2N = 64, the search 
did not converge even with A = 10~^. 

V. SUMMARY 

We have generalized the "Newton descent" variational method to determination of in- 
variant m-tori in general d-dimensional dynamical systems, and provided numerical evi- 
dence that the method converges in a large domain of existence of invariant tori, up to 
their breakups. In case of maps and flows with invariant tori such as standard maps, the 
approach offers an alternative method for determining critical thresholds. While in princi- 
ple the method is applicable to flows or maps in arbitrary dimension, computation can be 
expensive for invariant objects larger than 1- and 2-tori. We have utilized the smoothness of 



0.02 

-0.02 
-0.04 
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the fictitious time evolution to introduce acceleration schemes which improve the efficiency 
of the method. 

In our numerical work, we have implemented the method in the constant shift (4) 
parametrization, Fourier representation of an m-torus. Other discretizations could be bet- 
ter suited to specific applications. For instance, if an invariant torus is close to its critical 
threshold, representation of small fractal structures requires inclusion of slowly decaying high 
wavenumber Fourier modes, and a large number of Fourier modes is needed to obtain an ac- 
curate representation. Furthermore, the discretization points distribute very non-uniformly 
when close to criticality. In this limit, other non-constant shift parametrizations of the 
torus dynamics might be more appropriate. For example, our method is of modest accuracy 
compared to some of current studies of critical tori, in particular Haro and de la Llave [11] 
computation of critical tori to 100 digits precision. 

In periodic orbit searches we have found the Newton descent approach robust, and very 
useful for finding periodic orbits in high-dimensional phase-spaces where good guesses for 
multi-shooting Newton routines are hard to find [43, 44]. Examples worked out here suggest 
that the method is also a robust starting point for m-dimensional invariant tori searches. 
Once an approximate invariant torus is found by the Newton descent method, it can be used 
a starting guess for a high precision method, such as some of the currently used Newton's 
methods in Fourier space representations of invariant tori. 
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